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ABSTRACT 


The  essentially  non-oscillatory  (ENO)  finite  difference  scheme  is  applied  to  systems 
of  conservation  laws  ut  +  f(u)I  =  0  of  mixed  hyperbolic-elliptic  type.  A  flux  splitting 
f(u)  =  f+(u)  +  f~(u),  with  the  corresponding  Jacobi  matrices  having  real  and  pos¬ 

itive/negative  eigenvalues,  is  used.  The  hyperbolic  ENO  operator  is  applied  separately  on 
f+(u)x  and  on  f~(u)x.  The  scheme  is  numerically  tested  on  the  van  der  Waals  equation  in 
fluid  dynamics.  We  observe  convergence  with  good  resolution  to  weak  solutions  for  various 
Riemann  problems,  which  are  then  numerically  checked  to  be  admissible  as  the  viscosity- 
capillarity  limits.  We  also  observe  the  interesting  phenomena  of  the  shrinking  of  elliptic 
regions  if  they  are  present  in  the  initial  conditions. 
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1  Introduction 


The  system  of  conservation  laws 


iit  +  f(u),  =  0 

11(1, 0)  =  u°(x) 


(1.1) 


is  hyperbolic  if  the  Jacobi  matrix  has  real  eigenvalues  and  a  complete  set  of  eigen¬ 

vectors.  In  recent  years  there  have  been  a  lot  of  activity  in  designing  stable  and  accurate 
numerical  methods  for  solving  systems  of  hyperbolic  conservation  laws.  The  ENO  (essen¬ 
tially  non-oscillatory)  high  order  finite  difference  method,  [5],  [6],  [12],  [13],  is  one  of  the 
successful  approaches.  The  philosophy  of  ENO  schemes  is  to  use  upwinding  and  adaptive 
stencils,  based  on  the  local  “wind”  direction  (the  sign  of  the  relevant  eigenvalue)  and  the 
local  smoothness,  in  each  of  the  local  characteristic  fields.  ENO  schemes  can  resolve  the 
interactions  of  shocks  and  waves  and  complicated  wave  structures  quite  well,  according  to 
the  numerical  examples  for  the  hyperbolic  Euler  equations  of  a  polytropic  gas  in  [6],  [13]. 
ENO  schemes  are  formally  high  order  accurate,  measured  by  local  truncation  errors.  In  some 
cases  the  actual  order  of  accuracy  may  degenerate  due  to  the  frequent  switching  of  stencils  in 
the  linearly  unstable  regions  [9],  [14].  There  is  a  simple  modification  of  ENO  schemes,  with 
no  additional  computational  cost,  [14],  which  overcomes  this  accuracy  degeneracy  problem. 
However,  our  experience  seems  to  show  that  the  original  ENO  scheme  has  very  good  resolu¬ 
tion  for  nonlinear  systems  [13],  [14],  suggesting  that  the  modification  may  not  be  necessary 
in  such  cases. 

If  the  Jacobi  matrix  in  (1-1)  has  complex  eigenvalues,  the  system  becomes  elliptic. 
The  initial  value  problem  (1.1)  may  not  be  well  posed  in  the  elliptic  regions.  However,  in 
applications,  (1.1)  may  still  be  of  mixed  hyperbolic-elliptic  type.  Examples  include  equations 
in  fluid  dynamics  [17],  elasticity  [8]  and  the  partial  differential  equations  related  to  Lorenz 
systems  [7],  just  to  name  a  few.  In  this  paper  we  use  the  van  der  Waals  equation  in  fluid 
dynamics 


vt  +  p(w)x  =  0,  wt  -  vx  -  0 

u(i,  0)  =  u°(i),  w(x,  0)  =  w°(x) 


(1.2) 


with 


p(w)  = 


RT 


(1.3) 


w  —  b  w2 

where  R,T,a,  b  are  all  positive  constants,  for  our  numerical  examples.  See,  for  example,  [17] 
for  details.  Equation  (1.2)  corresponds  to  (1.1)  with  u  =  ( v,w)T ,  f(u)  =  (p(w),  -v)T .  The 
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two  eigenvalues  of  the  Jacobi  matrix  are  ±^/ —p'( w).  For  an  ideal  gas,  the  pressure  p{w ) 
is  a  decreasing  function  of  w,  resulting  in  a  hyperbolic  system  (1.2).  However,  during  the 
cu-existence  of  gas  and  liquid,  p'(w )  may  become  positive  within  an  interval,  as  in  the  case 
(1.3)  with  suitable  parameters  (Figure  1),  making  the  system  (1.2)  elliptic  in  this  region.  The 
mathematical  ill-posedness  of  the  system  in  the  elliptic  region  reveals  the  physical  fact  that 
the  state  in  the  elliptic  region  is  not  stable,  and  it  typically  evolves  into  “phase  transitions” , 
i.e.,  jumps  across  the  elliptic  regions  in  the  weak  solution.  As  in  the  hyperbolic  case,  there  can 
be  more  than  one  weak  solution,  and  there  are  efforts  in  the  literature,  e.g.,  [15], [16],  [8],  [2], 
[3],  [10],  [11],  to  establish  admissibility  criteria  with  the  goal  to  single  out  one  “physically 
relevant”  weak  solution.  A  weak  solution  u  =  ( v,w )  of  (1.2)  is  called  admissible  as  the 
viscosity-capillarity  limits  in  [15]  if  it  is  the  bounded  a.e.  limit  of  ue  =  (ve,wf),  satisfying 

'  Vt  +  P(w‘)x  =  evxx  -  e2Awxxx, 

<  w£t-vcx-0  (1.4) 

ve(x,  0)  =  v°(x),  we(x,  0)  =  w°(x) 

as  e  — >  0+.  Here  A  is  a  positive  constant,  usually  between  0  and  1.  The  analysis  usually 
starts  with  the  Riemann  problem 


[v(x,  0),  w(x,  0)) 


(vL,wL),  x<0 

(vr,  WR),  X  >  0 


(1.5) 


whose  solutions  contain  most  of  the  rich  features  (shocks,  phase  boundaries,  rarefaction 
waves,  etc.)  for  more  general  initial  conditions.  It  is  easier  to  consider  each  jump  (shock, 
phase  transition)  separately,  with  a  jump  called  locally  admissible  if  it  is  the  limit  of  the 
travelling  wave  solutions  of  (1.4).  A  weak  solution  is  then  called  locally  admissible  if  each 
of  its  jumps  is  locally  admissible.  Admissible  or  locally  admissible  solutions  of  the  Riemann 
problems,  e.g.  [15], [10], [3],  typically  contain  phase  boundaries  (discontinuities  jumping  across 
the  elliptic  regions)  but  do  not  achieve  values  inside  the  elliptic  region.  More  general  initial 
data  is  more  difficult  to  analyse.  It  is  interesting  to  know  what  happens  if  the  initial  condition 
contains  elliptic  regions.  ^From  a  computational  point  of  view,  if  a  shock  capturing  method, 
such  as  the  ENO  method  in  [6],  [13]  is  to  be  used,  discontinuities  are  typically  spread  out  in 
two  or  three  points,  hence  points  can  sit  in  elliptic  regions  even  if  the  exact  solution  jumps 
across  it. 


One  of  the  main  ingredients  of  ENO  schemes,  and  of  many  other  non-oscillatory  schemes 
such  as  TVD  (total- variation-diminishing)  schemes  [4],  is  the  approximation  in  each  of  the 
local  characteristic  fields.  If  the  system  becomes  elliptic,  local  characteristic  decomposition 
is  no  longer  available.  In  Section  2  we  propose  a  flux  splitting  f(u)  =  f+(u)  +  f~(u),  with 
the  corresponding  Jacobi  matrices  having  real  and  positive/negative  eigenvalues.  This 
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is  similar  to  the  flux  splitting  used  for  hyperbolic  systems,  for  example  the  Lax-Friedrichs 
splitting  and  the  van  Leer  splitting  [19].  The  hyperbolic  ENO  operator  is  then  applied 
separately  on  f+(u)I  and  on  f“(u)I,  using  the  component  version  in  the  elliptic  regions, 
since  the  characteristic  decompositions  of  have  no  physical  meanings  in  such  regions. 

This  idea  is  described  in  the  next  section.  In  Section  3  we  numerically  test  the  scheme  on 
the  van  der  Waals  equation  (1.2)-(1.3).  We  observe  convergence  with  good  resolution  to 
weak  solutions  for  various  Riemann  problems.  We  then  numerically  check  that  these  weak 
solutions  are  admissible  as  the  viscosity-capillarity  limits  by  computing  (1.4)  with  a  sequence 
of  decreasing  e.  We  also  compute  the  solution  with  smooth  initial  conditions,  and  observe 
the  interesting  phenomena  of  the  shrinking  of  elliptic  regions  if  they  are  present  in  the  initial 
conditions. 


2  The  Numerical  Scheme 

We  start  with  a  flux  splitting 


f(u)  =  f+(u)  +  f~(u)  (2.1) 

with  the  requirement  that  the  Jacobi  matrix  has  only  real  and  positive  eigenvalues, 

and  likewise  that  the  Jacobi  matrix  drg^  has  only  real  and  negative  eigenvalues.  If  the 
system  is  hyperbolic,  the  simplest  way  to  achieve  such  splitting  is  due  to  Lax-Friedrichs 


f±(u)  =  ^(f(u)  ±  au),  a  =  max  |  A,(u)  |  (2.2) 

where  Aj(u)  are  the  eigenvalues  of  the  Jacobi  matrix  For  special  classes  of  hyperbolic 

systems,  such  as  the  Euler  equations  of  a  polytropic  gas,  more  sophiscated  splittings  with 
better  physical  meanings  are  available,  e.g.,  van  Leer’s  splitting  [19].  For  an  elliptic  system, 
the  simple  splitting  (2.2)  no  longer  works.  In  fact,  any  splitting  with  the  Jacobi  matrices 
dra^  and  -fg^  commuting  with  each  other,  as  is  the  case  in  (2.2),  will  probably  fail, 
because  commuting  matrices  with  distinct  eigenvalues  can  be  simultaneously  diagonalized, 
hence  the  eigenvalues  of  their  sum  are  simply  the  sum  of  their  corresponding  eigenvalues. 
However,  a  splitting  similar  to  (2.2)  with  the  scalar  a  replaced  by  a  diagonal  matrix: 


( 


f±(u)  =  ^(f(u)  ±<*u),  a  = 


a2 


/ 


(2.3) 
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can  usually  yield  a  required  splitting.  For  example,  the  flux  f(u)  in  the  van  der  Waals 
equation  (1.2)  can  be  split  successfully  using  (2.3)  with: 

,  =  a2  +  M  (2.4) 

and 

M  =  a  max  y/max(0,p'(u;)),  a  >  2  (2.5) 

The  idea  is  to  make  an  ansatz  g(u)  =  M(v,  0)T,  then  try  to  find  the  smallest  possible  M 
such  that  the  Jacobi  matrices  d^u|^S(u))  both  have  real  and  distinct  eigenvalues.  This  leads 
to  M  given  by  (2.5).  Once  this  is  done,  it  is  easy  to  use  the  Lax-Friedrichs  idea,  i.e.,  to  add 
and  substract  au  with  a  suitable  a,  to  accomplish  the  splitting. 

Equipped  with  the  splitting  (2.1),  one  can  then  apply  any  successful  hyperbolic  approxi¬ 
mation  techniques  separately  to  f+(u)  and  f~(u).  The  only  exception  is  that  characteristic 
decompositions  should  not  be  performed  in  elliptic  regions,  since  the  characteristic  directions 
of  and  — do  not  have  any  physical  meaning.  In  this  paper  we  apply  the  ENO 

techniques  developped  in  [12],  [13]  to  f+(u)  and  f~(u). 

We  summarize  the  algorithm  briefly  in  the  following.  More  details  can  be  found  in  [12], 
[13]. 

(1)  The  spatial  operator 


a2  =  max  max  0, 


'  M2  —  4  p'(w)  —  M 


-  *W)x  =  -  f  (u)x 


is  approximated  by  a  conservative  flux  difference 


(2.6) 


L(“),  =  r+(u), +  L  (“), 

where  j  is  the  grid  index  (the  grids  are  located  at  Xj  —  j  Ax),  and  the  resulting  ODE 


Qt  -  L(u);  (2-8) 

is  discretized  in  the  time  variable  t  by  a  class  of  TVD  (total-variation-diminishing)  Runge- 
Kutta  type  high  order  methods  introduced  in  [12].  For  example,  the  third  order  case  is 


’  UW  =  u(°)  +  AfL(uW) 
u(2)  =  |u(°)  -I-  +  lAtL(u(^) 

uO)  =  lu(°)  +  |u(2)  +  |AfL(u^2)) 

,  u(°)  =  Un,  Un+1  = 


(2.9) 
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In  the  following  we  will  discuss  the  spatial  operator  L(u)^  in  (2.7)  and  suppress  the  time 
variable  t\ 

(2)  The  numerical  fluxes  f±±k  in  (2.7)  approximate  h±(xJ±x )  to  high  order,  where  the 
functions  h±(x)  are  defined  implicitly  by 


(2.10) 


see  [13]  for  the  derivation,  r-th  order  polynomial  interpolation  is  used,  based  on  Newton 
differences.  The  approximation  to  h±(x)  is  carried  out  component  by  component.  In  the 
following  we  will  use  plain  letters  to  represent  one  component  of  the  corresponding  vector  in 
bold  letters,  and  we  will  suppress  the  superscripts  ±.  As  is  pointed  out  in  [13],  there  is  no 
need  to  construct  the  function  h(x)  or  its  Newton  differences  explicitly:  one  can  simply  use 
the  difference  tables  of  /( u(x)).  If  the  undivided  differences  of  /(u(x))  are  computed  by 


flh  0]  =  /(t*i) 

=  f[j  +  1.*  -  1]  -  1]» 


k  =  l,...,r 


(2.11) 


where  r  is  the  order  of  the  interpolation  polynomial,  then  the  numerical  flux  /  +i  is  obtained 

by 


/,+ i  =  J2  -  h  *)/(»,  *] 


(2.12) 


where  i  is  the  left-most  point  in  the  stencil  used  to  approximate  /J+ 1 ,  and  C7(s,  k)  is  defined 

by 


I  »+fc  »+fc 


(2.13) 


The  small  matrix  C  is  independent  of  the  flux,  hence  can  be  computed  only  once  and 
then  stored.  The  reader  can  easily  check  that  (2. 12)-(2. 13)  gives  the  correct  interpolant 
approximation  to  the  function  h(x). 

What  distinguishes  ENO  from  other  finite  difference  methods  is  the  adaptive  stencil 
idea.  This  idea  is  realized  through  the  choice  of  i,  the  left-most  point  in  the  stencil  used 
to  approximate  /J+x .  We  start  with  i  =  j  for  computing  ft+L  or  i  =  j  +  1  for  computing 
f~+i-  This  is  due  to  upwinding,  since  the  eigenvalues  of  are  all  positive,  hence  the 

information  for  f+(u)  propagates  to  the  right,  likewise  the  information  for  f~(u )  propagates 
to  the  right.  The  stencil  is  then  expanded  point  by  point,  according  to  the  principle  of 
choosing  the  smaller  in  absolute  value  of  the  two  relevant  differences: 
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if  ( abs(f[i ,  k\).gt.abs(f[i  -  1,  k]))  i  =  i  —  1 


(2.14) 


for  k  =  1, r, 

A  factor  can  be  introduced  in  (2.14)  to  bias  towards  the  choice  of  a  linearly  stable  centered 
stencil  [14].  This  modification  can  enhance  the  accuracy  in  some  cases  but  may  cause  some 
over- compression  effects  [14].  Since  for  nonlinear  systems  the  original  ENO  scheme  works 
well  numerically  [6],  [13],  we  do  not  use  this  modification  in  this  paper. 

The  scheme  described  above  is  uniformly  high  order  in  space  and  time,  measured  by  local 
truncation  error  analysis. 

Remark  2.1  Since  the  schemes  described  above  are  conservative,  any  converged  solution 
will  be  a  weak  solution  of  (1.1).  It  is  more  difficult  to  show  that  the  limit  solutions  are 
admissible.  If  we  take  r,  the  order  of  the  interpolating  polynomial,  to  be  zero,  and  use  the 
splitting  (2.1)-(2.2)  for  a  hyperbolic  system,  we  recover  the  classical  Lax- Friedrichs  scheme. 
It  is  well  known  that  the  Lax-Friedrichs  scheme  can  be  rewritten  as  a  centered  scheme 
plus  a  dissipation  term  approximately  equal  to  |aAsu„.  If  we  use  the  splitting  of  the 
form  (2.3),  we  can  still  think  about  it  as  a  centered  scheme  plus  a  dissipation  term  of  the 
form  |aAiuIx.  It  is  then  reasonable,  cf.  [16],  [2],  to  expect  that  the  scheme  converges  to 
the  admissible  weak  solutions.  Ample  numerical  tests  should  be  performed  to  assess  the 
convergence  and  admissability  for  higher  order  schemes.  The  numerical  examples  in  Section 
3  are  preliminary  results  in  this  direction. 

Remark  2.2  If  some  fractional  step  method  (e.g.  Strang  [18])  is  used  on  the  splitting 
(2.1),  we  end  up  with  a  scheme  of  the  form  (see  (2.7)): 

u"+1  -  (I  AfL+)(I+AtL-)un  (2.15) 

(for  the  next  time  step  the  two  operators  may  reverse  order).  For  a  linear  or  nonlinear 
problem  with  smooth  solutions,  it  is  very  easy  to  choose  stable  operators  (I  +  AtL±),  due 
to  the  hyperbolicity  of  f±(u)  in  (2.1).  However,  this  does  not  necessarily  mean  that  the 
scheme  (2.15)  is  stable,  since  (I  +  AtL^1)  may  not  commute  with  each  other  and  may  not  be 
simultaneously  diagonalizable.  If  the  operators  satisfy  the  more  restrictive  condition 

||I  +  AfL+||  <  1  +0{At)  ||I  +  A<L"||  <  l  +  0(Ai)  (2.16) 

for  any  consistent  norm,  the  fractional  step  scheme  (2.15)  will  be  stable. 
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3  Numerical  Examples 

We  use  the  van  der  Waals  equation  (1.2)-(1.3)  with  RT  =  1,  a  =  0.9  and  b  =  0.25.  The 
graph  of  the  corresponding  p(w)  is  in  Figure  1.  The  system  is  elliptic  for  ct  <  w  <  0, 
where  a  =  0.574912  and  (3  =  1.036251.  The  so-called  Maxwell  line  BE  in  Figure  1,  where 
the  two  shaded  areas  are  equal,  intersects  the  curve  of  p(w)  at  w  =  m  =  0.494273  and 
w  =  M  —  1.405065.  The  horizontal  lines  AD  and  CF  in  Figure  1  yield  7  =  0.483100  and 
6  =  1.918618. 

We  use  the  third  order,  i.e.  (2.9)  and  (2.12)  with  r  =  3,  ENO  scheme,  described  in 
Section  2.  If  the  computational  cell  is  contained  completely  inside  one  of  the  hyperbolic 
regions  w  <  a  or  w  >  (3,  we  use  characteristic  decompositions  (the  ENO-LF  algorithm 
described  in  [12],  [13]).  Otherwise  the  component  by  component  approximation  described 
in  Section  2  is  used.  The  splitting  used  is  (2.3)-(2.5)  with  a  =  2.2.  The  time  step  A t  is 
restricted  by  a  CFL  number  0.6,  i.e.  At  <  0.6(p(~^-)  +  p(~^^-))Ax  where  p(A)  is  the 
spectral  radius  of  A. 

All  the  computations  are  performed  by  using  a  sequence  of  refined  meshes  to  verify 
convergence,  although  we  typically  only  show  the  graphs  for  one  or  two  fixed  meshes. 

We  first  compute  several  Riemann  problems  (1.5): 

(1)  ( vl,wl )  =  (l,m),  ( vr,wr )  =  (1,  M)  where  m  and  M  are  the  Maxwell  values  de¬ 
fined  above.  This  initial  condition  satisfies  the  Rankine-Hugoniot  condition  for  a  stationary 
jump.  Physical  principles  (Maxwell  equal  area  rule)  and  many  admissibility  criteria  (e.g. 
the  viscosity-capillarity  criterion  in  [15],  see  (1.4))  indicate  that  this  is  an  admissible  jump. 
Our  numerical  result  shows  a  stable,  sharp  jump  for  this  case,  Figure  2.  We  remark  that 
here  and  in  what  follows,  the  numerical  solution  usually  has  one  or  two  transition  points 
in  the  elliptic  region  for  the  phase  jump.  Apparently  they  do  not  cause  any  trouble  to  the 
computation. 

(2)  (vi,wi)  =  (1,0.54),  {vr,wr)  =  (1,0.54).  This  initial  condition  also  satisfies  the 
Rankine-Hugoniot  condition  for  a  stationary  jump,  but  physical  principles  and  many  admis¬ 
sibility  criteria  (e.g.  the  viscosity-capillarity  criterion  in  [15],  see  also  [10])  indicate  that  this 
is  not  an  admissible  jump.  Our  numerical  result  shows  the  evolution  of  this  jump  into  a 
more  complicated  structure  of  jumps,  Figure  3a,  apparently  due  to  the  inherent  numerical 
viscosity  of  the  scheme  (see  Remark  2.1).  The  solution  exibits  oscillatory  behaviors  near  the 
phase  boundary,  Figure  3a.  This  is  unpleasant  but  not  surprising  since  we  used  component 
by  component  approximations  in  cells  involving  elliptic  regions,  hence  during  the  process  of 
one  wave  splitting  into  two  or  more  waves,  one  or  more  of  them  being  hyperbolic,  oscillations 
occur  as  a  failure  of  recognition  of  the  corresponding  characteristic  fields.  Similar  oscillations 
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also  appear  for  hyperbolic  systems  if  a  component  by  component  approximation  is  used  (see, 
e.g.  [1]).  The  oscillations  become  smaller  and  more  confined  when  the  number  of  grids  is 
increased  (Figure  3b),  indicating  the  convergence  of  the  scheme  even  with  these  oscillations. 

The  background  of  Figure  3a  is  computed  by  the  same  scheme  with  2000  points.  It  agrees 
with  the  result  with  4000  points  hence  can  be  considered  as  a  converged  solution.  In  order 
to  check  whether  this  weak  solution  is  admissible  as  a  viscosity- capillarity  limit,  we  plot 
in  Figure  3c  the  numerical  solutions  of  (1.4),  for  A  =  with  e  =  0.1,0.01,0.001  and  the 
solution  of  our  scheme  for  (1.2).  The  solutions  to  (1.4)  are  computed  by  the  standard  fourth 
order  centered  scheme  with  the  classical  fourth  order  Runge-Kutta  time  discretization.  We 
verify  adequate  resolution  for  the  solution  of  (1.4)  for  each  fixed  e  by  repeatedly  refining  the 
mesh  until  the  solutions  do  not  change  to  visual  inspection  (the  largest  number  of  grid  points 
used  is  8000).  Clearly  we  can  see  the  convergence  of  the  solutions  of  (1.4)  to  our  solution 
when  e  — >  0+  in  Figure  3c. 

(3)  —  (1,0.45),  (vr,Wr)  =  (2,1.5).  This  case  is  somewhat  easier  to  compute  than 

the  previous  case,  since  the  initial  condition  is  not  a  steady  nonadmissible  weak  solution. 
Figure  4a  shows  the  result  with  200  grid  points  on  a  background  of  a  converged  solution 
with  2000  grid  points.  Figure  4b  shows  the  convergence  as  e  — ►  0+  of  the  solutions  of  the 
viscosity-capillarity  equation  (1.4)  to  our  solution. 

We  then  compute  the  solutions  for  smooth  initial  conditions: 

(1)  (u°(i),u;0(x))  =  (1  —  0.5cos(x),  1  +  0.5sm(x)).  This  initial  condition  crosses  the 
elliptic  regions  (Figure  5).  The  solution  gradually  evolves  into  piecewise  smooth  solutions 
contained  entirely  inside  one  of  the  two  hyperbolic  regions  w  <  a  and  w  >P,  connected  by 
jumps  over  the  elliptic  regions  (phase  transitions).  This  seems  to  agree  with  the  physical 
intuition. 

(2)  (u0(x),io°(x))  =  (1  —  0.5cos(x),  0.8  +  0.2sin(x)).  This  initial  condition  is  entirely 
contained  in  the  elliptic  region.  However,  similar  to  the  previous  case,  the  solution  gradually 
evolves  into  piecewise  smooth  solutions  separated  by  phase  transitions,  Figure  6.  Notice  that 
during  the  evolution  oscillations  are  generated  inside  the  elliptic  regions,  presumably  due  to 
the  inherent  instability  of  the  equation  in  those  regions.  These  oscillations  fade  out  once  the 
solution  evolves  into  the  hyperbolic  regions. 


4  Concluding  Remarks 

f 

Numerical  methods  for  solving  systems  of  conservation  laws  of  mixed  hyperbolic-elliptic  / 
type  are  investigated,  through  a  flux  splitting  to  write  the  physical  elliptic  flux  as  a  sum 
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of  two  hyperbolic  fluxes  with  positive/negative  eigenvalues,  then  to  apply  the  essentially 
non-oscillatory  (ENO)  high  order  finite  difference  methods  on  each  of  them.  The  method, 
in  the  simplest  first  order  case,  is  equivalent  to  adding  a  numerical  dissipation  term  with 
a  diagonal  dissipation  matrix.  The  numerical  results  on  the  van  der  Waals  equation  of  gas 
dynamics  indicate  that  the  method  can  resolve  phase  boundaries  well  and  can  be  used  as  a 
tool  to  study  the  evolution  of  elliptic  regions.  More  numerical  tests  on  different  mixed  type 
equations  constitute  current  research.  ■'  /  f  /  _____ _ „ 


i 
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Figure  1:  p(w )  = 

0.55 
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Figure  2:  Maxwell  solution,  (vL, wL)  =  (l,m),  ( vr,wr )  =  (1,M),  t=25.  200  points  (plus) 
and  the  exact  solution  (solid  line).  Here  and  in  what  follows,  the  region  between  the  two 
dashed  lines  is  elliptic. 


CX  n3 


Figure  4:  (vL,wL)  =  (1,0.45),  ( vR,wR )  -  (2,1.5),  t=1.5.  4a:  200  points  (plus)  and  2000 
points  (solid  line);  4b:  centered  solutions  of  (1.4)  with  e  =  0.1,  0.01,0.001,  and  ENO  solution 
for  (1.2). 

Figure  4a 


Figure  4b 
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